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We present an implementation of an efficient algorithm for the calculation of the spectrum of 
one-dimensional quantum systems with periodic boundary conditions. This algorithm is based 
on a matrix product representation for quantum states (MPS), and a similar representation for 
Hamiltonians and other operators (MPO). It is significantly more efficient for systems of about 100 
sites and more than for small quantum systems. We apply the formalism to calculate the ground 
state and first excited state of a spin-1 Heisenberg ring and deduce the size of the Haldane gap. The 
results are compared to previous high-precision DMRG calculations. Furthermore, we study spin-1 
systems with a biquadratic nearest-neighbor interaction and show first results of an application to 
a mesoscopic Hubbard ring of spinless Fermions which carries a persistent current. 



I. Introduction 

It was recognized early on that density matrix 
renormalization group (DMRG) simulations of one- 
dimensional (ID) quantum systems require significantly 
more numerical resources for periodic boundary con- 
ditions (PBC) than for open boundary conditions 
(OBC) % Verstraete, Porras, and Cirac (VPC) ad- 
dressed this issue, and they proposed an algorithm in 
terms of matrix product states (MPS), which scales sig- 
nificantly better with the matrix size to of the MPS than 
standard DMRG with PBC. However, intermediate steps 
of this algorithm require matrices of size to x to 2 , and 
computer time and memory necessary to determine the 
improved representation still scales with to 5 as compared 
to to 3 for OBC. 

This issue was addressed by Pippan, White, and Ev- 
ertz (PWE) @, who recognized that for sufficiently large 
systems a much more efficient implementation is possible 
using a singular value decomposition (SVD) of products 
of certain transfer matrices. In order to calculate such 
products with sufficient accuracy only rather few singu- 
lar values must be kept. 

The usefulness of the improved algorithm was demon- 
strated in Ref. by a calculation of the ground state of 
the spin-1 Heisenberg Hamiltonian. The authors showed 
that accurate results for the ground state energy are ob- 
tained by a comparison with highly accurate standard 
DMRG calculations. As a result, it was concluded that 
for large enough systems one obtains an algorithm which 
scales similarly with to as calculations for systems with 
OBC. 

In the present paper we extend the PWE algorithm in 
two respects: First, we propose an implementation of this 
algorithm in terms of MPS and matrix product operators 
(MPO). To this end we define generalized transfer matri- 
ces, which are subjected to an SVD. This enables further 
gains in efficiency in certain situations. Second, we ex- 



tend the PWE framework and include the calculation of 
excited states of ID many body Hamiltonians. 

We apply this algorithm to a small selection of spin 
models (bilinear and biquadratic spin-1), as well as to a 
spinless Fcrmion model. In the course of these applica- 
tions it was found that in general the number of singular 
values one must keep depends on the matrix size to, i.e. 
the larger to the more singular values must be kept in 
order to produce high precision results. 

From the MPS representation it is straightforward 
to calculate correlation functions and other observables. 
Results of such calculations will be presented elsewhere. 



II. MPS-MPO formalism for PBC 

We first rewrite the algorithm proposed in Ref. Q in 
terms of MPS and MPO: The states of a ID quantum 
system of size N (e.g. a spin system) are approximated 
in terms of a matrix product state (MPS), 
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Here, the Oj represent the local degrees of freedom at 
site j, and each represents a matrix of size to x to, 
where m is called bond dimension. In the algorithm to be 
described the elements of these matrices are variational 
parameters to be adjusted using a suitable optimization 
procedure. The trace in Eq. {1} ensures periodic bound- 
ary conditions and includes a sum over all aj. 

Analogously, operators are written as matrix product 
operators (MPO) 



<U Wi,...,a N )(a[,...,a' N \, (2) 



and the trace includes a sum over all <jj and o'y 



gam, 



each W^ a , represents a matrix of size my/ x my/, i-e. 
each W is a tensor of order 4. It turns out that all op- 
erators of interest with short range interactions (e.g. the 
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Heisenberg Hamiltonian) can be written in terms of W 
tensors with small bond dimensions mw- The structure 
of the W tensors is determined by the specific model 
under investigation. Wc will provide the explicit MPO 
representation of the various operators later in in this 
paper. 

Matrix elements of MPO in such states, 

(0|O|V>> = Tr E$ (A, B) ■ . . . ■ E$ (A, B), (3) 

can be expressed in terms of the (generalized) transfer 
matrices 

E® (A, B)=Y, W ™' ® )* ® B ® ■ ( 4 ) 

a, a' 

The matrices A and B characterize the states \4>) and \ijj), 
respectively. The Kronecker product <& in Eq.|4]obviously 
produces transfer matrices of size m 2 mw x m 2 mw- For 
later use we also define the special transfer matrix 

eW(a,b) = Y,&°AaWt®bW. (5) 

One advantage of the MPO formalism used here over 
the formalism employed by VPC and PWE is the fact 
that it takes care of the structure of the effective Hamil- 
tonian to be determined automatically (as encoded in the 
MPO), while the effective Hamiltonian in the VPC for- 
mulation depends structurally on the Hamiltonian of the 
model under consideration. 

In order to find the ground state of the many body 
system one solves a standard variational problem using 
the matrix elements of the MPS as variational parame- 
ters. The optimization of the variational parameters of 
the MPS is implemented as a local update step, which is 
repeated until convergence is achieved @. In the MPO 
formalism such a local update step amounts to the solu- 
tion of a generalized eigenvalue problem 

H%<pV\ = e b1jvjp[i] (6) 

in terms of the effective Hamiltonian H e g and the effec- 
tive normalization matrix N e g given by 

m w / - — ■ — \ 

nj& = Y, w ki®( H n- H n > ( 7 ) 

kl v / lk 

ivj = E ® hff-N^j . (8) 

The energy of the state is obtained from and this 
value will converge to the ground state energy eventu- 
ally. In fact, we stop the iterative update procedure, if 
this quantity does not change any more with respect to 
defined convergence criteria. 

The updated MPS is obtained from tpU] = M( CTi j/ n by a 
suitable partitioning of the vector into a tensor. The tilde 
in (|7|) indicates the operation X(ij\uiji} = -2t(ii'),rM') for 
each m 2 x m 2 submatrix of the bracketed quantities. As 



a consequence of this transposition the effective Hamil- 
tonian and the normalization matrix are assured to be 
Hermitian matrices and standard methods for the solu- 
tion of generalized eigenvalue problems can be applied. 
(For open boundary conditions the normalization matrix 
is unity and only a standard eigenvalue problem needs to 
be solved.) 

The matrices H^} , and H [ £ , N$ are the prod- 
ucts of transfer matrices from all sites to the left and 
to the right of the site j, where the MPS is updated. 
The H matrices are obtained from generalized transfer 
matrices as defined in Eq. (fj|, while the N matrices are 
formed from the transfer matrices defined in Eq. (|5|), in 
both cases setting A = B = M with M the MPS to be 
determined. 

In the algorithm proposed by VPC one sweeps back 
and forth over the entire lattice several times updating 
the MPS at each site until convergence of the energy 
e^l is achieved. Initially, one starts from a randomly 
selected MPS. After each update step the updated matrix 
is regauged in order to keep the algorithm stable. The 
standard regauging procedure, which assures the relation 

J2bWbW = 1 (9) 

a 

after each update step is described in more detail in 
Refs. and 0. 

Similarly, excited states will be constructed itcratively 
by finding the lowest state in the space orthogonal to 
the space spanned by the states already found. Wc will 
denote the matrices of these MPS by $^' fe where k enu- 
merates these states (k = for the ground state, k = 1 
for the first excited state, etc.). It was pointed out in 
Ref. Q that this construction can also be implemented 
itcratively as an update step by locally projecting to the 
orthogonal subspace. Here, we need to determine the 
local projection operator with the property 

pblybl = o v k (10) 

with 

Y k ] = °r ■ [ [ ] <S> [ 1) and Y^-Y^} = if k ± m. (11) 

Here, the spin and m indices of the 9z. k matrices are 
suitably combined to form a vector. For simplicity, we 
will use the same symbol $b1 for these vectors (see the 
analogous definition of 0^1 above). 

The matrices O^} and O^ are products of transfer ma- 
trices as defined in Eq. (JSJ) from all sites to the left and to 
the right of the site j, respectively, and setting B = M 
and A = $ fc with M the (excited) MPS to be determined. 
The update procedure for these matrices is implemented 
as a generalized eigenvalue problem (see Eq. ©) for the 
projected effective Hamiltonian PbljJ eff pblt ) anc [ nor _ 
malization matrices pb'l jv eff pb']t_ The- (local) projection 
operator Pbl will be constructed according to Eq. ([TO")) 
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by finding a set of vectors orthogonal to the calculated 
y fe . A standard numerical orthogonalization routine is 
employed for this purpose. 

III. Efficient implementation 

In order to implement the local update steps just de- 
scribed one needs to calculate various products of trans- 
fer matrices. These are standard matrix products, which, 
however, depending on the bond dimension of the MPS 
and MPO, they may be numerically expensive. Naively, 
a multiplication of two transfer matrices Q requires 
0(m mf^) operations, which may be reduced in view of 
the structure of the transfer matrices to 0(m 5 m^,). In 
analogy to the proposal by PWE we will now describe a 
procedure to reduce this operational count further. This 
reduction occurs due to the structure of the W tensors 
and, in particular, for products of transfer matrices with 
many factors, i.e. long products. Here (unlike Ref. 0]) we 
consider products of transfer matrices in terms of MPSs 
and MPOs, 

E$(A,B)-...-E®(A,B)= J2 *k ( 12 ) 

fc=i 

As was pointed out by PWE the sum over k may be 
cut at rather low values, which for the generalized trans- 
fer matrices has two reasons: First, the rank mj of the 
transfer matrices is in many practical situations lower 
than mwfn 2 ■ This reduces the upper limit of the sum to 
ms- E.g. as is indicated below, the rank of the transfer 
matrices for the Ising or Heisenberg models is 2m 2 and 
not 3m 2 or 5m 2 , respectively, as expected naively. This 
reduction of the summation limit is exact and does not 
depend on the product length. 

However, for long products, the upper limit may be 
reduced to very low values due to the fact that only very 
few singular values in the expansion Eq. (fT2"j) are sig- 
nificantly different from 0. For ground state calculations 
of chains with about 100 sites and m = 10 one needs to 
consider only about 20 singular values. This is demon- 
strated for the Heisenberg model in Fig. [TJ This figure 
corresponds to Fig. 1 of Ref. [§J and shows rather sim- 
ilar results for the Nl. Here, we also plot the singular 
values of Hl , and we see that only a few more singular 
values than for Nl are needed. (Beyond a certain limit 
the singular values are set to an irrelevant small constant 
by our computer implementation.) 

In order to utilize this feature for the local update al- 
gorithm described in the previous section one needs to 
implement the algorithm in such a way, that only suffi- 
ciently long products of transfer matrices occur during 
the update process. Therefore, one cannot use the stan- 
dard sweeping procedure since 'short' products of trans- 
fer matrices occur at the turning points of the sweeps. 
Following PWE we implement the algorithm as a circu- 
lar update procedure. The ring of sites is separated into 



three sections as shown in Fig. [2J and the update pro- 
cess occurs always in the 'active' section. The algorithm 
is then implemented in 3 basic steps: 

1. (Initialization step) Start from some initial ran- 
domly created matrix product state \ip) as defined 
in Eq. [TJ The bond dimension of all matrices 
(j = l,...,N) is m. Partition the set of matri- 
ces into three sections as shown in Fig. [3J 
Initialize section 3 with a singular value decompo- 
sition (SVD) of the products of generalized transfer 
matrices defined in Eqs. (U]) and and store this 
SVD in the tensors hi and nl, respectively. Ini- 
tialize section 2 with an SVD of the products of 
transfer matrices and store this SVD in the tensors 
hr and nr. 

2. (Update step) Goto section 1. Initialize each site of 
section 1 with the appropriate product of transfer 
matrices moving counter-clockwise starting from 
the product corresponding to section 2. Then up- 
date and regauge the MPS in section 1 moving 
clock-wise using the previously calculated products 
of transfer matrices. Updating means solving the 
generalized eigenvalue problem described above for 
each site. (One immediately obtains an SVD of the 
products of transfer matrices inside the updated 
section by multiplication to the SVD of the previ- 
ous site, i.e. one does not need to calculated an 
SVD at each update step. This is an important 
advantage of the algorithm using MPS and MPO.) 

Finally copy the tensors nl and hi on the tensors 
nr and hr and calculate the SVD of the product of 
transfer matrices of section 1 with the just updated 
MPS matrices and store this SVD in the tensors nl 
and hi. 

3. Goto section 2 and do analogous calculations as de- 
scribed for section 1 above. Continue with further 
steps moving clockwise to the neighboring section 
until convergence is achieved. 

An important prerequisite for the implementation of 
the algorithm is an efficient SVD. This has been described 
in Ref. [H, but we have a few remarks: Let M be a prod- 
uct of transfer matrices. Then, according to the proce- 
dure outlined in Ref. @ one has to form products of these 
matrices M with some matrices x and y' of size p x m 2 , 
e.g. y = xM and z = My' T . In order to do this effi- 
ciently one must not calculate the matrix M explicitly, 
but rather multiply each transfer matrix in M recursively 
to x or y' starting from one or the other end of the se- 
quence of factors in M. Then the multiplication of M to 
the matrices x or y' can be done in 0(Npm 3 ), where N 
is the number of transfer matrices in M. 

Similar steps as outlined above for ground state cal- 
culations are required for the determination of excited 
states, i.e. for each excited state we use the same algo- 
rithm searching for the optimal MPS in the space orthog- 
onal to the states already found. We have implemented 
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FIG. 1: Distribution of the singular values of products of transfer matrices for m=10, Nl (left) and Hl (right), with 33 terms, which is 
the minimum number of terms in our calculation for N = 100 sites on a ring with homogenous nearest-neighbor Heisenberg interactions. 



Construction of the MPO for periodic boundary con- 
ditions is not difficult, 




FIG. 2: Circular algorithm for a ring with N sites: the ring is 
partitioned into three sections. Updating only happens in one of 
them, so that we always deal products of transfer matrices with 
minimum length N/3. For further discussion see the main text. 



the described algorithm within a few pages of Mathemat- 
ica code. 



IV. Matrix product operators 

In order to apply the algorithm developed above to 
specific problems we must define the relevant degrees of 
freedom, the size of the local Hilbert space, and the in- 
teraction in terms of a suitable MPO. Once this MPO is 
defined, the implementation of the algorithm takes care 
of the details of the calculation. 

The simplest examples to be considered are spin mod- 
els, e.g. the spin-S' unisotropic Heisenberg Hamiltonian 
in an external magnetic field B, 

N 

n = jJ2s?®s? +1 + s?®sy +1 

i=l 

N 

+AS*®S* +1 -BY / S?, (13) 



with the exchange interaction J, and the unisotropy A. 
In the following we will set J = 1. The Hamiltonian is 



written in terms of the spin operators Si 
l 



jcr,-, and for 



2 L 

spin- ^ the a matrices correspond to the standard Pauli 



matrices. Periodic boundary conditions correspond to 
setting N + 1\->1. 
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for i 



,N 



with e a unit matrix. The local single-body Hilbert space 
has dimension 2S+1, and the bond dimension is dw = 5. 
However, the rank of the transfer matrices which deter- 
mines the cost of the calculation is not 5m 2 as expected 
naively but only 2m 2 . The first matrix has a different 
structure as the other matrices in order to fulfill the re- 
quired boundary conditions. 

For a bilinear-biquadratic spin-S" ring with the Hamil- 
tonian 



H = aSi ® sv 



b(S t <g> 



^+i) 2 



(15) 



one easily finds an explicit MPO representation with a 
bond dimension dw = 14. Here, again, the rank of the 
transfer matrices is not 14m 2 but only 2m 2 , which re- 
duces calculational cost significantly. 

Calculation of matrix elements for observables (e.g. 
the magnetization or correlation functions) is straight 
forward in the MPS representation either with an MPO 
representation of the operators or without. Also for these 
calculations one may take advantage of the fact that such 
calculations are just products of transfer matrices (see 
(Eq.[3])) and use the expansion (TT2"|) for long products. In 
the present paper we will use this feature for the calcu- 
lation of the variance of the Hamiltonian as is discussed 
in the next chapter. 
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V. Applications 

In order to test the implementation of the proposed 
algorithm we start out with calculations of the isotropic 
Heisenberg model also studied in Ref. Q. Of course, it 
is easily possible to calculate energy spectra for small 
systems (10-50 sites) using our implementation, and we 
have calculated up to 30 excited states for such systems. 
However, then one must take into account most or all of 
the singular values in the expansion of the transfer matri- 
ces. In order to take advantage of a significant reduction 
of the number of singular values, the system size should 
be about 100 sites or more, and we present results for 
systems with 100 sites in this paper. 

In order to run such calculations three important pa- 
rameters, which determine the precision of the results 
must be set: The bond dimension m, the number of sin- 
gular values to be included in the expansion of the various 
transfer matrices p and p', and the number of update 
steps N u , where p is the number of singular values re- 
tained in the expansion of the Nx matrices, and p' those 
of the Hx matrices. 

Of course, a large to is desirable, however, the algo- 
rithm scales at least with p'm 3 N, so in practice we are 
presently limited to about to = 50. We shall demonstrate 
below, that the number of singular values to be taken into 
account increases with m, and one must be careful not to 
take too few terms in the expansion Eq. (|T2|) . Unfortu- 
nately, convergence of the update process is rather slow 
close to the minimum of the energy. Therefore, for high 
precision results we need more and more update steps, 
and usually we choose their number dynamically by ob- 
serving the change of the calculated energy within one 
sector. If this change (averaged over the whole section) 
is below a certain limit, we stop the update process. 

One purpose of the present calculations is to gain ex- 
perience which parameter setting for m,p,p', and N u is 
required in order to find e.g. the Haldane gap in a spin-1 
ring with sufficient precision. The gap is obtained from 
a difference of two large energies of similar value, so the 
two energies must be calculated with rather good preci- 
sion. (Let us note parenthetically that the m required 
in the present algorithm is significantly smaller than the 
corresponding quantity in standard DMRG calculations.) 

In Fig. [3] we show the distribution of singular values 
of the transfer matrices obtained at the end of a calcu- 
lation with m = 30 for the isotropic spin-1 Heisenberg 
model, i.e. the calculations shown in Fig. [T] and Fig. [3] 
only differ in the choice for to. From a comparison of 
these results one concludes that if one increases m one 
also needs to increase the number of singular values to 
be taken into account. Our experience shows that the 
necessary increase is quite significant depending on the 
m one wants to use for a particular calculation. This fact 
was not mentioned in Ref. Q. However, after this paper 
was nearly completed, we became aware that a similar 
observation was made in Ref. [f| for the standard PWE 
algorithm without MPO. 



The MPS-MPO formalism employed here allows to 
straightforwardly test how well the calculated MPS ap- 
proximates an eigenstate of the Hamiltonian. To this end 
one calculates the variance 

AH = (H 2 ) - (H) 2 (16) 

which should be zero for an eigenstate. Since from the al- 
gorithm we obtain an explicit representation of the state, 
we can, at least in principle, easily evaluate this quantity, 
if we find a suitable MPO representation of H 2 . The bond 
dimension of H 2 is m?^ , but its rank is often significantly 
lower, which is used to significantly reduce the cost for 
the calculation of H 2 using the expansion Eq. (fT!?)) . 

The results obtained so far for the isotropic Heisenberg 
model are summarized in Table HI The ground state en- 
ergy is in good agreement with that reported in Ref. [|| . 
In addition we show results for the first excited state from 
which wc determine the Haldane gap, which agrees with 
the infinite system DMRG calculations of Ref. [l| to two 
significant digits. Haldane @ conjectured on the basis 
of a field theoretical study that gencrically integer spin 
chains are gapped, while half-integer spin chains are gap- 
less in the thermodynamic limit. For specific examples 
(spin-i and spin-|) we can confirm this numerically with 
our calculations. 

For the ground state, we judge the precision of the ob- 
tained results by a comparison to a high precision result 
obtained by DMRG as quoted in Ref. Q, and assume 
that this value is numerically exact for the Heisenberg 
ring with 100 sites. In fact, this result is quite close to 
the infinite system value obtained in Ref. [1(. 

A second interesting test of the implementation of the 
proposed algorithm is the biquadratic chain Eq. (|15p 
(with a = and b = —1) investigated in detail by 
S0rensen and Young [?J using a mapping of the bi- 
quadratic spin-1 ring to the XXZ spin-i system, which 
can be solved using Bethe Ansatz techniques. In Table [TXI 
we present some preliminary results for this system using 
our technique, which are compared to the high-precision 
Bethe Ansatz results of Ref. [7| . In the thermodynamic 
limit one expects a doubly degenerate ground state and 
a small gap to the next excited state. Of course, for fi- 
nite systems the degeneracy is lifted. This system is an 
interesting testing ground for our numerical techniques 
as there are extremely precise results available for sys- 
tems up to 1000 spins. Only for such large systems one 
expects to be close to the thermodynamic limit. 

The results indicate good agreement with the Bethe 
Ansatz results, however, for high precision one needs 
large to and for m = 30 one needs about 30-60 singu- 
lar values to be taken into account. Convergence of the 
energies at a particular to, depending on the precision 
required, may be slow. Therefore, we recommend to cal- 
culate first with a few different m in order to see the to 
dependence before one iterates with the chosen m to high 
precision. How well the calculated MPS approximates an 
eigenstate is measured by a calculation of AH. 

As a last example we apply the proposed algorithm 
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TABLE I: Ground state energy Eo, first excited state energy E\, and Haldane gap Ei — Eo for an isotropic spin-1 Heisenberg 
ring of N = 100 sites. A/E is the relative difference between our calculated result and the value calculated by DMRG given in 
Ref. @]. 



m 


Eo/N 


A/E 


Ei/N (3) 


Ei — Eo 


10 


-1.40122726344 


1.83 10" 4 


-1.39621210860 


0.50153 


20 


-1.40145874749 


1.47 10" 5 


-1.39730198769 


0.41566 


30 


-1.40148324293 


5.83 10" 7 


-1.39736419879 


0.41192 


40 


-1.40148390219 


9.73 10" 8 


-1.39737237500 


0.41115 


DMRG [3] 
DMRG (infinite) [1] 


-1.4014840386(5) 
-1.40148403897 




-1.39737901875 


0.41050 



Length=33 



Length=33 



30 
k 



FIG. 3: Distribution of the singular values of products of transfer matrices for m=30, Nl (left) and Hl (right), with 33 terms, which is 
the minimum number of terms in our calculation for N = 100 sites on a ring with homogenous nearest-neighbor Heisenberg interactions. 



to a Hubbard model of spinlcss Fcrmions, and in par- 
ticular to a mesoscopic ring filled with electrons pierced 
by a magnetic field, such that persistent currents can be 
observed. The Hamiltonian of this system is given by 



N N 

H = -tJ2 (Ac e+1 e-^ /N +h.c. S j+U^2n i n t+1 +Vn 1 



resentations of the single-body operators read 



£=1 



(17) 



Here </> is the magnetic flux piercing the ring, U the 
nearest-neighbor Coulomb interaction and V the local 
interaction of an impurity at site 1. Here, ct and c are 
Fcrmion creation and destruction operators, and n the 
density operator. The hopping energy t will be set to 1, 
and periodicity requires to set N + 1 n- 1 . 

More details about this Hamiltonian and its physics 
may be found in Refs. and the references therein. The 
Hamiltonian is U(l) symmetric, and the particle number 
is a good quantum number to label the states. Due to 
the impurity, the model is not homogeneous: it is one 
advantage of our MPS implementation that it can han- 
dle inhomogeneous problems, since it does not assume 
translational invariance of system. 

Since we are considering spinless Fermions, the local 
single-body Hamiltonian describes a two-level system, 
which is analogous to a spin-^ system. The matrix rep- 
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(18) 



Together with the 2x2 unit matrix these matrices (like 
the Pauli matrices) form a complete set. 

One then obtains the following MPO representation 
for this problem, 
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in terms of the matrices defined in Eq. (|18j) . the param- 
eters of the Hamiltonian, and a chemical potential [i to 
be discussed below. The minus signs in the first row of 
arise due to the anti-commutativity of the creation 
and destruction operators on different sites. 
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TABLE II: Ground state energy Eo, first excited state energy Ei and gap E\ — Eq for a biquadratic spin-1 Heisenberg ring of 
N = 100 sites (a = 0, b = —1 in Eq. (|15[) ). A/E is the relative difference between the Bethe Ansatz results and the numerical 
values obtained, (p = 30, p' = 60), AH the variance of the Hamiltonian 



m 


Eo/N ( A/E ) 


AH 


Ei/N ( A/E ) 


AH 


Ei — Eq 


10 


-2.794 020 092 (1.04 10" a ) 


1.08 


-2.793 830 121 (1.05 10" 3 ) 


1.10 


0.018 997 


20 


-2.795 792 099 (4.07 10" 4 ) 


0.44 


-2.795 632 899 (4.12 10~ 4 ) 


0.44 


0.016 077 


30 


-2.796 790 186 (5.03 10" 5 ) 


0.03 


-2.796 675 842 (3.95 1(T 5 ) 


0.28 


0.011 452 


Bethe Ansatz [7] 


-2.796 930 734 




-2.796 786 305 




0.014 442 



In order to study persistent currents one needs to cal- 
culate the ground state energy as a function of the mag- 
netic flux and then calculate the persistent current j us- 
ing the Hellmann-Feynman theorem, j = —dE ((f>)/d(f>. 

Since experiments are usually made for systems with 
fixed particle number, it would be necessary to develop 
the algorithm in such way that it respects the U(l) sym- 
metry of the Hamiltonian. At this stage our implemen- 
tation does not respect this symmetry. Of course, it is 
always possible to shift to the ground state of the sector 
with the desired particle number using an appropriate 
chemical potential fx. However, this chemical potential is 
usually not known, and one would need to use an itera- 
tion process to find that chemical potential such that the 
resulting state contains the desired number of particles. 
Only for half-filled systems, it is known that the required 
chemical potential to find the ground state equals the in- 
teraction U. We therefore concentrate here on half-filled 
systems, and shift the spectrum accordingly. 

First results are shown in Table IIIII for a ring with 
N = 128 sites. In order to be able to calculate persistent 
currents using Hellmann-Feynman theorem one must be 
able to precisely distinguish the ground state energies 
for different 0, which requires rather high-precision cal- 
culations. The energy determined for the ground state 
agrees with the result given in Ref. We also calcu- 
late the energy of the next higher/lower state and the 
number of particles n it contains. Clearly, the ground 
state is half-filled, while the next higher/lower state con- 
tains one particle less. At cf> = one finds a degenerate 
ground state in the half-filled sector. (Here, our proce- 
dure to calculate 'excited' states, may yield even a lower 
lying state, since within the spectrum there exist states 
below the ground state of the half filled sector.) For fu- 



ture calculations an implementation respecting the U(l) 
symmetry is desirable. 

TABLE III: Energy of the half-filled ground state E and 
energy of the next higher/lower state Ei of a spinless Fermion 
ring described by the Hamiltonian Eq. (|17|l for N — 128, 
m = 30, U = 1, and V = 0. 



<t> 


Eo 


n 


Ei 


n 





-63.98647233 


64 


-63.98581164 


64 


tt/2 


-64.00411240 


64 


-64.94361781 


63 


7T 


-64.01004832 


64 


-64.94770847 


63 



VI. Summary 

In this paper we suggest a new version of an efficient 
MPS algorithm for one dimensional systems with peri- 
odic boundary conditions. The present version unlike 
the original proposal Q uses an MPO representation. 
We also extend the algorithm for the calculation of ex- 
cited states. We report about first results obtained with 
this algorithm, and investigate the necessary parameter 
settings in order to obtain high precision results for sys- 
tems with 100 sites. The advantage of the algorithm is 
that one obtains an explicit representation of the many- 
body quantum state, which can then be used to calculate 
observables such as correlation functions. We will report 
about such calculations in a forthcoming publication. 
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